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A NON-SYMMETRIC COUPLING OF THE FINITE VOLUME METHOD 
AND THE BOUNDARY ELEMENT METHOD 

CHRISTOPH ERATH, GUNTHER OF, AND FRANCISCO-JAVIER SAYAS 

Abstract. As model problem we consider the prototype for flow and transport of a con¬ 
centration in porous media in an interior domain and couple it with a diffusion process in 
the corresponding unbounded exterior domain. To solve the problem we develop a new 
non-symmetric coupling between the vertex-centered finite volume and boundary element 
method. This discretization provides naturally conservation of local fluxes and with an 
upwind option also stability in the convection dominated case. We aim to provide a first 
rigorous analysis of the system for different model parameters; stability, convergence, and 
a priori estimates. This includes the use of an implicit stabilization, known from the finite 
element and boundary element method coupling. Some numerical experiments conclude the 
work and confirm the theoretical results. 

Keywords, finite volume method, boundary element method, non-symmetric coupling, 
convection dominated, existence and uniqueness, convergence, a priori estimate 

Mathematics subject classification. 65N08, 65N38, 65N12, 65N15 


1. Model problem and introduction 

Throughout this work, let fl C R d , d = 2, 3, be a bounded domain with connected polygonal 
Lipschitz boundary T and O e = is the corresponding unbounded exterior domain. We 

consider the same model problem as in |Eral21 IEral3a] : find u and u e such that 


AVu + b u) + cu 

= / 

in 0, 



(la) 

-A u e 

= 0 

in f2 e , 



(lb) 

U e (x) 

= Coo log \x\ + 0{l/\x\) 

for x -A oo, 

d = 

2, 

(lc) 

U e (x) 

= 0(l/\x\) 

for x -A oo, 

d = 

3, 

(Id) 

u 

= u e + u 0 

on T, 



(le) 

(AVu — bu) ■ n 

du e 

— o + *0 

an 

on r n , 



(If) 

(AVu) • n 

du e 

— a + *0 

an 

on T out , 



(lg) 


where A is a symmetric diffusion matrix, b is a possibly dominating velocity field, c is a 
reaction function, / is a source term, and Coo is an unknown constant. The coefficients are 
allowed to be variable. The coupling boundary T = dVt = dfl e is divided in an inflow and 
outflow part, namely T m := {x G T | b(x) • n(x) < 0} and T out := {x G T | b(x) • n(x) > 0}, 
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respectively, where n is the normal vector on T pointing outward with respect to G. We 
allow prescribed jumps Uo and to on T. The radiation condition for the two dimensional case, 
which will be complemented with the additional hypothesis that the diameter of is less 
than one, guarantees that our problem has a unique solution. Other radiation conditions 
are also possible, but some lead to restrictions on the data. Changing from one to the other 
is a relatively simple exercise adding sources. See [McLOOl 1CS85J for more information on 
radiation conditions. 

The model problem in the interior domain hi is the prototype for flow and transport of a 
concentration in porous media. Usually, boundary values such as Dirichlet and/or Neumann 
boundary conditions are needed to solve the problem. These problems are often convection 
dominated and the conservation law, e.g., local conservation of fluxes, should also be pre¬ 
served for a numerical approximation of the solution. Therefore, a finite volume method 
(FVM) is often the method of choice since they provide an easy option to stabilize the 
convection term and they natural preserve conservation of numerical fluxes due to their for¬ 
mulation. However, if the domain is unbounded one would have to truncate the domain. 
The above formulation solves also another issue, i.e., if we do not know any boundary con¬ 
ditions, we assume a diffusion process in the corresponding (unbounded) exterior domain 
fl e , which “replaces” the boundary values. The method of choice for unbounded domains 
is the boundary element method (BEM) which reduces the discretization to the boundary 
and therefore avoids the truncation of Vt e . Therefore, we consider an FVM-BEM coupling as 
in |Eral01 IEral2l IEral3aj . To the best of the authors knowledge, these works are the first 
theoretical justifications of a FVM-BEM coupling, where a three held coupling approach is 
used with either the vertex-centered (finite volume element method, box method) FVM or 
the cell-centered FVM. 

In this work we analyze and verify a non-symmetric FVM-BEM coupling with the vertex- 
centered FVM, in the following only named FVM. The main motivation of using this is to 
get an easier coupling formulation and a smaller system of linear equations, which saves com¬ 
putational costs. The idea of a non-symmetric coupling approach goes back to [ JN80 1 , fBJ79 j. 
This coupling formulation applied for a finite element method (FEM)-BEM discretization is 
also known as Johnson-Nedelec coupling. However, the analysis in this early works relied 
on specific choices of the discretization spaces or on the compactness of a certain integral 
operator, which was in fact a restriction to a smooth boundary. In particular, a rigorous 
mathematical analysis for Lipschitz domains was not known. Recently, the work in | |Say09| 
provided a first analysis, which overcame these restrictions. Meanwhile, several extensions 
and simplifications are possible, such that a SIAM review paper |Sayl3| was published. 
Among these extensions there are results on the non-symmetric formulation for the poten¬ 
tial equation with variable coefficients |QS131 IStellj , non-linearities |AFF + 13l IFFKP15] , for 
elasticity |FFKP15~1 IStel3j . and for boundary value problems |GHS12l IOS14] , In addition, 
similar results have been reported on related coupling formulations |AFF + 131IGHS12] and the 
DG-BEM coupling [ HS15j . We want to mention that the counterpart to the non-symmetric 
coupling is the so called the symmetric coupling first introduces in |Gos87] . However, sym¬ 
metry is referred to a diffusion-diffusion transmission problem, i.e., the whole system is 
symmetric. We stress that this would be destroyed if one applies convection in the interior 
domain. 
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There exist a couple of papers, which analyze the vertex-centered FVM, e.g., [ BR871 lHac89j 
to mention only the very first works. It is well known that for pure diffusion with piecewise 
constant diffusion coefficient on a primal mesh the standard FEM and the FVM bilinear form 
are exactly the same. Thus the schemes differ basically only on the right-hand side. However, 
for all other diffusion problems |Cai91] and a possible convection field and a reaction term 
the systems are different. Contrary to standard FEM, FVM still provides local flux conser¬ 
vation due to its formulation and provides an easy upwind stability option for convection 
dominated problems. The standard analysis approach makes use of a comparison between 
the FEM and FVM bilinear form [BR871 IHac891 ICai911 IELL021 ICha02] . For our FVM-BEM 
coupling we may apply similar techniques for the FVM part. Note that contrary to a classical 
FEM-BEM coupling we do not have a classical Galerkin orthogonality property due to the 
FVM formulation based on the conservation law. Thus the analysis differs significantly to 
an FEM-BEM analysis. However, we use the equivalent formulation of a stabilized continu¬ 
ous coupling formulation, extended here for the convection-diffusion-reaction problem in H, 
and compute an ellipticity constant. Based on the continuous stabilization we introduce a 
stabilization for the FVM-BEM coupling. This is needed for pure diffusion models and for 
convection-diffusion-reaction problems, where the energy norm reduced to a semi-norm. We 
stress that the stabilization is only needed for theoretical purposes since the formulation is 
equivalent to the standard system. We aim to provide a discrete ellipticity estimate, conver¬ 
gence, and a priori estimates for the FVM-BEM coupling. Our new analysis technique gives 
us a recipe for the coupling of BEM with a non-Galerkin method like FVM. Furthermore, 
this work improves the results in jEralOl lEra!2] for a three field FVM-BEM coupling, where 
we had to assume a little bit more regularity on the unknown exterior conormal solution and 
some constraints on the convection and reaction terms for some special model problem con¬ 
figurations. However, as for the non-symmetric FEM-BEM coupling we have a theoretical 
constraint on the eigenvalues of A, which is not needed in (EralOl IEral2j . 

Throughout, we denote by L m (-) and H m (-), m > 0 the standard Lebesgue and Sobolev 
spaces equipped with the usual norms || • ||z, 2 (-) and || • ||ff m (-), respectively. For w C O, (-, ■)„ is 
the L 2 scalar product. The space /7 m_1 / 2 (r) is the space of all traces of functions from H m (Q) 
and the duality between H m (T) and is given by the extended L 2 -scalar product 

(-, -)r- The space H} oc (Vt ) := {u : H —y M | v\k £ i7 1 (RT), for all K C Hopen and bounded} 
collects functions with local H 1 behavior. Furthermore, the Sobolev space W 1)OQ contains 
exactly the Lipschitz continuous functions. If it is clear from the context, we do not use a 
notational difference for functions in a domain and its traces. To simplify the presentation 
we equip the space T~L := i7 1 (H) x H~ l / 2 ( T) with the norm 

IMk := IMIffRn) + 11 ^11#- 1/2 (r) 


for v = (v, VO £ 'H- 

With this notation we can specify the model data as: the diffusion matrix A : — y M. dxd has 
entries in W 1 ’ 00 (f2), is bounded, symmetric and uniformly positive definite, i.e., there exist 
positive constants (7a, i and Ca ,2 with (7a,i|v| 2 < v 7 A(x)v < (7a,2|v| 2 for all v e and 
almost every x £ O. We will also admit coefficients A that are T-piecewise constant, where 
T denotes the triangulation of H introduced in subsection 13.11 satisfying identical symmetry 
and uniform positive definiteness assumptions. Note that the best constant Ca.i equals the 
infimum over x £ fi of the minimum eigenvalue of A(x), which we will denote A m i n (A). 
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Furthermore, b G W 1,00 (fi) d and c G L°°(Q) satisfy 

j(x) := — divb(x) + c(x), 'j(x) > 0 for almost every x G 

with the function 7 G L°°(Q). We stress that our analysis holds for constant b and c = 0 as 
well. Finally, we choose the right-hand side / G L 2 (Q), u 0 G if 1//2 (r), and t 0 G H~ l / 2 (T). In 
the two dimensional case we additionally assume diam(Q) < 1 to ensure H~ l / 2 (T) ellipticity 
of the single layer operator defined below. 

Then our model problem reads in a weak sense: find u G and u e G HLJ£ L) such 

that ©-© hold. 

The model problem (JTJ) admits a unique solution for both, the two and three dimensional 
case [Eral2j . 


Remark 1. To replace the radiation condition (fTc|) by u e (x) = 0{l/\x\) for |x| —> 00 in two 
dimensions one would have to assume the the scaling condition 

(du e /dn, l) r = 0 

to guarantee solvability. As opposed for the purely diffusive case, this condition cannot be 
easily transformed into a condition on the data. 


The content of this paper is organized as follows. Section 2 gives a short summary on integral 
equations and the weak formulation of our model problem based on the non-symmetric 
approach. We show an ellipticity estimate through an equivalent stabilized weak formulation 
and state the ellipticity constant explicitly. In section 3 we introduce the non-symmetric 
FVM-BEM coupling to solve our model problem. Section 4 proves stability, convergence, 
and an a priori result for our coupling. Numerical experiments, found in section 5, confirm 
the theoretical results. Some conclusions complete to work. 


2 . Integral equation and weak coupling formulation 

The representation formula for the exterior Laplace equation (llb|) with the radiation condi¬ 
tion (fTc|) - (ITdl) and <f{x) = -$-u e (x) |r, iGK reads 


u e (x) = - G(x - y)4>{y) ds y + 


/ — 

Jr Jr & n y 

with the fundamental solution for the Laplace operator 


G(x - y)u e (y)\ r ds v 


(3) 


G{z) : = 


_ *r log|z| 
1 1 
47T Id 


for z G M 2 \{0}, 
for z G M 3 \{0}. 


(4) 


From (J3]) we obtain (taking traces) the boundary integral equation on T 

u e |r = (1/2 + JC)u e \r - Vf. 

The single layer operator V and the double layer operator /C are given, for smooth enough 
input, by 

(Vf)(x)= I ip(y)G(x - y) ds y (0)(x) = f 9(y)-^-G(x - y) ds y iGT, 
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where r\. y is a normal vector with respect to y. The integral equation (J5J) holds on T except 
on corners and edges. We recall | ICos 88 i Theorem 1] that these operators can be extended 
to bounded operators 

VeL(if a_ 1 / 2 (r);if a+ 1 / 2 (r)), ice L(H S+1/2 {T)-H S+1/2 (T)), se[-|,±]. 

It is also well-known that V is symmetric and H~ l / 2 (T) elliptic, since we additionally assume 
diam(Q) < 1 in the two dimensional case, which can always be achieved by scaling. The 
expressions 

II • llv := <Vv>r, || • ||v-i := 

dehne norms in i7 _ 1 / 2 (T) and H l / 2 (Y), respectively. These norms are equivalent to the usual 
ones. 

We consider a weak form of the model problem (jT]) in terms of boundary integral operators. 
For that we use the non-symmetric approach, i.e, calculate the weak formulation of the 
interior problem and replace the interior conormal derivative by the exterior 0 := du e jd n|r 
and the corresponding jump relations to, OB (|lg|). Second, we take the weak form of (TXT) 
and replace the exterior trace u e \r by the interior trace W|r and the jump u 0 , (ITcj) . Then the 
coupling reads: find u E if 1 (Q), 0 E i/" 1 / 2 (T) such that 

A(u, v) - (0, v) r = (/, v) n + (t 0 , v) T , (5a) 

(0, (1/2 - K)u)t + (0, V0) r = (0, (1/2 - /C)tt 0 )r (5b) 

for all v E // 1 (h2), 0 E H~ l ^ 2 (T). The bilinear form in (I5a| is given by 

A[u, v ) := (AVm — b u, Vr)n + (cu, v)n + (b • n u, v)r°^- 

Lemma 2. The bilinear form A is coercive and continuous on i/ 1 (h2) x H 1 ( fl), i.e., for all 
v,w E if 1 (J2) and y(x) from assumption (J2]) there holds 

{ ^'U,ilMI# 1 (fi) for l{ x ) > 0 almost everywhere in hi, 

^AilMlirqn) for 'y(x) >0 on ui C 0, |u;| > 0 , 7 ( 0 ;) = 0 elsewhere , ( 6 ) 

for y(x) = 0 almost everywhere in Cl, 

\A(w , v )I < CU.olHItfMnilMIffHn)- (7) 

Here, the constants C AA = min{A min (A), inf xen 7 ( 0 ;)} > 0, C' A X = A min (A) > 0 and C AA > 
0, depend on the data A, b and c. The constant C A1 = min{A min (A), C( i y(x), uj, O)} > 0 
depends additionally on the constant C( , y(x),uj,Tl) > 0, which is not known but depends on 
7(x) > 0 in uj, uj, and hi. 

Proof. There holds 


/ b • n v 2 ds > - 

I p out 2 


bn v 2 ds = - / div(bu 2 ) dx = -((div b)u, v)q + (by, ’Vv)q. 


If Adivb(a;) + c(x) > 'j(x) > 0 of assumption 
follows that 


is positive almost everywhere in hi, it 


A(v,v) > (AVv, Vv)q + -((divb)u,u) n + (cv,v) n > C AA \\ 
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If 7(2;) > 0 holds on a set co C f2 of positive measure but 7(2;) = 0 on Q\oj, we can use 
a compactness argument (or the Deny-Lions theorem) to prove coercivity of A in Z/ 1 (0). 
Then the coercivity constant l is not known. When 7(2;) = 0 almost everywhere in hi, we 
only obtain coercivity of A with respect to the H 1 seminorm and the constant C\ x . Using 
simple arguments, the continuity bound ([7]) can be easily proved with 

Ca ,2 — 2 max{|| A|| i <x>(^)dxd, 1111(ri)j ||c||l°°(Q)} + Crl|b • nlU^r™ 4 ), 

where Cr is the norm of the the trace operator i/ 1 (h2) — y L 2 (T out ). □ 

For convenience the system (l5al) - (15bD can be written in the product space T~L = x 

H-V 2 (r) as follows: we introduce the bilinear form B : U x U —> R 

B{(u : 0); (v, ip)) := A(u, v) - {<j>, v) r + (0, (1/2 - K)u) T + (0, V0) r , (8) 

and the linear functional 

F((v, ip)) := (/, v) n + (to, «)r + {1>, (1/2 - JC)u 0 ) r . (9) 

Then (l5a ]) - (|5b|) is equivalent to: hnd u E H such that 

£>(u; v) = F(v) for all v e %. (10) 

With integration by parts we calculate 

£?(v; v) = (AVr, Vu)n + ((| divb + c)v, v)q — (b ■ nv, v) r ™ + (b • nr, v) r out 
~ {ip, v) r + {ip, (1/2 - K)v ) r + (0, V0) r , 


and thus we see 

£>(( 1 , 0 ); ( 1 , 0 )) = / (1 div b + c) + / |b-n|. 

Jn Jr 

Thus if 1 div b + c = 0 in U and b • n = 0 on T (in particular, when b = (0, 0) T and c = 0), it 
follows that £>((1, 0); (1, 0)) = 0. This lack of coercivity will be remedied using an equivalent 
variational problem for the sake of analysis. 

Therefore, we define the linear operator 


P((v, ip)) ■= (1, ( 1/2 - K)v + Vip )r = / (( 1/2 - K)v + V0) 


/ 3 : = 


and introduce a parameter (3 depending on 7 ( 2 :) of assumption (J 2 ]); 

1 if 'y(x) = 0 almost everywhere in U, 

0 else. 

Then the /3-dependent perturbations of the bilinear form B( u, v) is 

B{ u;v) := £?(u, v) + /3P(u) P(v), 

and of the linear map F(v) 

F(v) := F(v) + /3(1, (1/2 - /C)u 0 )rP(v). 


(ii) 


( 12 ) 

(13) 


Thus a stabilized variational formulation is given by: hnd u G P such that 

B (u; v) = F(v) for all v 6 Ji. 
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(14) 


Note that this type of stabilization has also been considered in [OS13] and [AFF +13] . We 
emphasize that this formulation is introduced purely for theoretical purposes, and the dis¬ 
cretization will be applied directly on (j5a]) - ([5b]) . 


Lemma 3. The variational formulation (HU and the stabilized version in (IFfl) are equivalent. 

in |AFF + 13 


Proof. The equivalence of formulations was stated in |AF F + 13 . Theorem 14] for a pure 
diffusion problem. The convection and reaction terms in the bilinear form A{-, •) do not 
affect the proof. We note that we will see a similar result for the FVM-BEM discretization 
in Lemma [12] □ 


The next theorem on the coercivity of the bilinear form B is an extended and improved 
version of the one stated in |OS131 Theorem 3.1] and [AFF + 13l Theorem 15] for a purely 
diffusive problem. We extend it by the convection and reaction terms in the bilinear form and 
present an improved ellipticity constant compared to |OS13l Theorem 3.1]. This is possible 
due to some modification of the proof inspired by |OS14j . Before we state the theorem, we 
recall an important contractivity result for the double layer operator [OS131 Lemma 2.1] 
with the contraction constant C/c from |SW01] : there exists C/c G [1/2,1) such that 


||(l /2 + IQvWI-, < C k (V~\1/2 + JC)v,v ) r . 

Furthermore, we define for /3 = 0 

{ inf j{x) for 'y(x) > 0 almost everywhere in f 2 , 

for y(x) > 0 on w C O, |u;| > 0 , 7 (x) = 0 elsewhere 


(15) 


(16) 


with y(x) from assumption (J 2 ]) and the unknown constant C(j{x), oj, fl) > 0 introduced in 
Lemma [2] 


Theorem 4. If A min (A) > C/c/ 4, then B is PL-elliptic. More precisely, for all v = (v, ip) G PL 
holds 


£(v;v) > Cstab [||Vu|| i2(n) 


The stability constant C s t a b reads 




(17) 


Cstab 


min jCbc, \ A min (A) + 1 - yj (A min ( A) - l ) 2 + C K | for p = 0, 


min |l, \ A min (A) + 1 - a/(A min (A) - l ) 2 + C K 
and depends on the model data A, b, c, and the contraction constant C £. 


for p — 1 


Remark 5. The right-hand side in lira defines an equivalent norm in PL. While this is 
obvious for fi = 0, a simple compactness argument (see |AFF + 13l Lemma 10 and (65)] for 
a similar argument) shows the equivalence for P — 1. Note that we only do not know the 
constant CAab explicitly in the second case of (fT6l) . 


Proof. The proof is in the spirit of previous publications [OSl3l ()S 11 , [A FF + T3] on the non- 
symmetric FEM-BEM coupling, but extended here for the different interior model problem. 
Therefore, we only present the key points. 
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An element v G i/ 1 (f2) can be decomposed as a sum v = vr + vq, where up is harmonic and 
no G Thus (Vur, Vw)n = 0 for all w G Hq(Q), which implies that 

l|Vn||| 2(n) = ||Vno||| 2(n) + ||Vn r ||| 2(n) = ||Vn 0 ||i 2(n) + (S^v,v) r , (18) 

where S mt := V _1 (l/2 + /C) denotes the Steklov-Poincare operator, i.e., the Dirichlet to 
Neumann map of the interior Laplace problem. The term ( S mt v , v)r will help to compensate 
possible negative contributions of the non-symmetric coupling to the total energy of the 
system. Let us first recall our choice of /? depending on r y(x) in (HIT) and the definition of 
Cbc in (USD- This allows us to write the coercivity estimate of Lemma [2] as 

A(v,v) > A min (A)|| Vn||^ 2(f2) + (1 - /3)C' bc ||^||| 2( n ) for all v G H\Q). 

Following [OS13j and using (ITS]) , we can easily estimate 

(tp, (1/2 + JC)v) T = (ViP, V~\l/2 + JC)v)r 

< ||V -1 (1/2 + K)v\\vUWv = ||(l/2 + /C)n|| v -i||^||v 

for all (v,ip) G TL. Therefore, for all v = ( v,ip) G TL, we can estimate 
B (v; v) = A(v, v ) + (ip, Vip ) r - (ip, (1/2 + JC)u) T + /3P(v) 2 

> A mta (A)||V«||i l(n) + (1 - /3)C bc ||t,||i s(!1) + /3P(v) 2 + ml - Cl /2 (Sm,v)l /2 |Mlv 

> A m m(A)||Vuo|| l 2 (0) T (1 — (3)Cbc\\v\\ 2 L2 ^ + /3P(v) 2 

(S mt v,v)]/ 2 \ ( A min (A) — \\/Ck\ f(S int v,v)\! 2 \ 

Mv ) K-WC .c 1 ) \ |H|v )’ 

where in the last inequality we have used the harmonic splitting (ITS]) . Since A min (A) > 0, 
the quadratic form in the right-hand side of the above estimate is positive definite if and 
only if 

Amin(A) — \\[Ck 

— \VCk 1 

Calculating the smallest eigenvalue of the matrix above, we can bound 

B (v; v) > C„» b (||VtJ 0 ||| a (!!) + (S mt v, v) T + (1 - /3)H»||| a(1! ) 

+ pP{y ) 2 + Wilt). 

which, using f[I5|) . is the estimate of the statement of the theorem. □ 

Remark 6. Note that this result also improves the estimate of [OS131 Theorem 3.1] for a 
pure diffusion problem in fh The smallest eigenvalue in CAab in the case j3 = 1 is observed 
to be sharp in the numerical experiments o/m®, contrary to the constant reported therein. 

Using the boundedness of A in (J7J) and mapping properties of the integral operators, it is 
easy to conclude that the bilinear form B defined in (fI2|) and the linear form F in (ITS]) are 
bounded. Thus we can conclude the unique solvability of (TTT|) . Due to the equivalence of the 
formulations in Lemma El the original variational formulation (I5a[) - (l5bl) is uniquely solvable. 


— Amin (A) — -jCf c > 0. 











Remark 7. Note that the equivalence of (TT4|) and ca shown in Lemma\3\and the ellipticity 
estimate era also hold true on the discrete level, if the constants are in the discretization 
space of In other words, a possible FEM-BEM coupling solution, as shown in 

Remark 1771 exists and is unique and the Cea Lemma applies. 


3. A NON-SYMMETRIC FVM-BEM COUPLING 

In this section we develop a FVM-BEM coupling discretization in the sense of a non- 
symmetric coupling approach. From now on we assume to £ F 2 (r). First, let us introduce 
the notation for the triangulation and some discrete function spaces. 

3.1. Triangulation. Throughout, T denotes a triangulation or primal mesh of 0, and A f 
and £ are the corresponding set of nodes and edges/faces, respectively. The elements T G T 
are non-degenerate triangles (2-D case) or tetrahedra (3-D case), and considered to be closed. 
For the Euclidean diameter of T G T we write hx sup xyeT \x — y\. Moreover, He denotes 
the length of an edge or Euclidean diameter of E G £. The triangulation is regular in 
the sense of Ciarlet |Cia78j . i.e., the ratio of the diameter hx of any element T G T to 
the diameter of its largest inscribed ball is bounded by a constant independent of hx, the so 
called shape-regularity constant. Additionally, we assume that the triangulation T is aligned 
with the discontinuities of the coefficients A, b, and c of the differential equation (if any), 
the data /, uq, and to. Throughout, if n appears in a boundary integral, it denotes the unit 
normal vector to the boundary pointing outward the domain. We denote by £x C £ the set 
of all edges/faces of T, i.e., £x '■= {E G £ | E C dT } and by £r := {E G £ \ E C Y } the set 
of all edges/faces on the boundary T. 


Dual mesh. We construct the dual mesh T* from the primal mesh T as follows. In two 
dimensions we connect the center of gravity of an element T G T with the midpoint of the 
edges E G £x\ see Figure Qfa), where the dashed lines are the new boxes, called control 
volumes. In three dimensions we connect the center of gravity of an element T G T with 
the centers of gravity of the four faces E G £x- Furthermore, each center of gravity of a face 
E G £x is connected by straight lines to the midpoints of its edges. The elements of this 
dual mesh T* are taken to be closed. Note that they are non-degenerate domains because 
of the non-degeneracy of the elements of the primal mesh. Given a vertex a* G Af from the 
primal mesh T (i — 1... ^A f), there exists a unique box containing a, . We thus number the 
elements of the dual mesh Vi E T*, following the numbering of vertices. 


Remark 8. In two dimensions, instead of starting the construction of the boxes in the center 
of gravity, we can use the center of the circle circumscribed to the element. Connecting 
these points with the midpoints of the edges we form the so called Voronoi or perpendicular 
bisector meshes, since the connection between to neighbor’s circumscribed circle points is 
perpendicular to the shared edge. Our analysis works with such meshes as well. 

Discrete function spaces. We define with S 1 (T) :={t)G C(f2) \v\ T affine for all TgT} 
the piecewise affine and globally continuous function space on T. The space V°(£r ) is 
the £r -piecewise constant function space. On the dual mesh T* we provide V°(fT*) 
|«G L 2 (Ct) | v\y constant V G T*}. With the aid of the characteristic function y* over the 
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Figure 1 . The construction of the dual mesh T* from the primal mesh T in 
two dimensions with the center of gravity point in the interior of the elements 
the dashed lines (gray boxes) are the new control volumes V) 

AT. In Figure (b) 


in Figure 
of T* and are associated with 


we see an example 
intersection r 17 = V\ D V 7 ^ 0 of two neighboring cells Vi, V 7 £ T*, where Ti 7 
is the union of two straight segments. For a 3 , a 4 £ A/", where both a 3 and a 4 
lie on T, t 34 = V 3 fl V 4 7^ 0 is only a single segment. 


volume Vi we write for £ V°(T*) 

* X' ^ * * 

v h= £^ v i Xi, 

XiGAf 

with real coefficients v t . Furthermore, we define the T*-piecewise constant interpolation 
operator 

T h : C(H) V°(T), T h v := v( ai )x*(x). (19) 

CLi GA/" 

Because of the construction of the dual mesh from the primal mesh and the definition of 
there hold the well known results: 


Lemma 9. Let T £ T and E £ St- For Vh £ 5 1 (T) there holds 


[ (v h - X* h v h ) ds = 0, 

J e 

(20) 

||Vh — 1*h V h II L 2 (T) < h T \ Vw/jl \l 2 (T): 

(21) 

IK — 2JK|| l 2 (e) < Ch]l 2 \\XJ Vh\\L 2 {T)i 

(22) 


where the constants C > 0 depend only on the shape regularity constant. 


Proof. The proofs are standard. Note that for (I2U|) we need the fact, that the dual mesh 
T* is constructed through the midpoint of an edge E £ S in the two dimensional case and 
the center of gravity point if E is a face in the three dimensional case. A proof of (I2T|) can 
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be found in [Era ID!, and (122] | follows from (12T| through the standard trace inequality. Note 
that the above statements are independent of the choice of the interior point in T G T for 
the T* construction. □ 


3.2. The discrete system. A classical finite volume discretization describes numerically 
a conservation law of the model problem, i.e., a quantity in a volume can only change due 
to the inflow and outflow flux balance through its boundary. More precisely, for our model 
problem we integrate (flajl over each dual control volume V G T* and apply the divergence 
theorem. If we use the transmission condition (fTfjW lg) we thus get a balance equation for 
the interior problem 


ldV\T 


(—A\7uh + buh)-nds+ / cuhdx 


'v 


(23) 


b • n Uhds 


I dvnv° ut 


4>hds= fdx+ t 0 ds 


'dvnr 


iv 


idvrr 


for all V G T*. Note that the discretization in the interior domain follows along the dual 
mesh T*. Here, Uh G 5 1 (T) and fih G V°(£r) approximate u and q 1, respectively. We can 
rewrite (1231) in terms of a variational formulation; 


Av(uh,Vh) — {(fth,ZhVh)r — (f,Zh v h)n + {to,Xfvh)r 


with the finite volume bilinear form Ay '■ 5 1 (T) x 5 1 (T) —> R given by 


A v (u h ,v h ) ■= ^2 v h(ai) 

ai&M 


'9Vi\r 
cuh dx + 


A Vuh + b Uh) ■ n ds 
b • n Uhds 


'V 


18Vinr out 


(24) 


Remark 10. Note that the trial and test spaces are different in practice. The test functions 
in the finite volume part are in which is realized by taking nodal values Vh^afi in (12^ 

and by interpolation Tfvh G V°(T*) for Vh G 5 1 (T). We have chosen the above definition to 
simplify the notation below. 


To complete the coupling formulation we choose as in the classical non-symmetric FEM- 
BEM formulation the BEM equation (j4]) and replace the continuous ansatz and test spaces 
by discrete subspaces. Finally, the discrete system reads: End Uh G 5 1 (T) and fih £ V°(£ r) 
such that 


A v (u h ,v h ) - {fi hl X* h v h ) r = (f,Th v h ) n + (t^l* h v h ) T , (25a) 

{4>h, (1/2 - JC)u h ) r + (ifh, Vfih) r = (i>h, (1/2 - JC)u 0 ) r (25b) 

for all v h G S\T), G V°(£ r ). 

As in the continuous case we write the system in a more compact way. We consider the 
product space Ti h := <S' (T) x V°(£r), the bilinear form By : TLh x Ti h —> R 

B v ((w h ,(j) h y, (v h ,fj h )) :=A v (w h ,v h ) - (fi h ,lfv h ) r 

+ {fi’h, (1/2 — IC)wh)r + {ifh, V(j>h)r, 


it 





and the linear functional Fy : Tih —> K 


Fy((vh, iph)) (f,Fh v h)n + (to,Fh v h)r + (Ad (1/2 — /C)w 0 )r- 
The f|25ap ~ fl25bp is equivalent to: find E FLh such that 

By(u h ] v fc ) = F v (v h ) for all v h e U h . 


(26) 

(27) 


3.3. Upwind scheme. In general it is a non trivial task to get a stable discrete solution 
for convection dominated problems. Finite volume schemes, however, allow an easy upwind 
stabilization [ IRST96] . If we want to apply an upwind scheme for the finite volume scheme, we 
replace b u h on the interior dual edges/faces V)\r in Ay (123]) by an upwinded approximation. 
Given Vi E T*, we consider the intersections with the neighboring cells r*j = Vi D Vj A 0 for 
Vj E T* . Note that in two dimensions T t] is the union of two straight segments or (when the 
associated vertices g, : , clj E AT lie on T) a single segment; see Figure Q [j/b)] In three dimensions 
Tij consists of one or two polygonal surfaces. We then compute the averages 


Aj : ~ 



F 

^3 

32 

i- 

A* := ,/ . / 

1 J Tij 

1 Aj 1 j Ti 


A ds, 


where ip points outward with respect to Vi , and the parameter 


Aj • ® (Aj | Fj | /1| Ay || oo), 

for a weight function $ : M —y [0,1], which is being applied to the Peclet number. Then we 
consider the value 


Uh,,i] • AijUhifli) + (1 Aj) " j h Aj) 


instead of Uh when restricted to t, 3 C dV t \ F. In this work we choose the upwind value 
defined by the classical (full) upwind scheme by 


<F(t) := (sign(f) + l)/2, 

i.e. A ij = 1 for Aj A 0 and A^ = 0 otherwise. A second choice will be 


*(*) == 


mm 


{2|t| _1 ,l}/2 for t < 0, 

1 — min {2|f| -1 , l}/2 for t > 0, 


(28) 


(29) 


where we can steer the amount of upwinding to reduce the excessive numerical diffusion. 
Whenever we apply an upwind scheme for the convection part, we replace the finite volume 
bilinear form Ay in the system (I25aj) - (125bp by 


Ay p (u h ,v h ) := v h( a i) 

ai&M 


-AVw/ l -nds+ / cuhdx 


' 8Vi\r 


'Vi 


/ b • n Uh,ij ds + b • n 

tfr. Jru ’ JdVinr™* 


Uh ds 


j£Afi JTi i 


(30) 


where denotes the index set of nodes in T of all neighbors of a* E AT. 
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4. Stability and convergence 


In this section we want to introduce a stabilized FVM-BEM coupling version of (12?)) for 
analysis purposes only. As in (fT4|) we use the “implicit theoretical” stabilization of [AFF +13] , 
Similar as above we define By : Ph x Ph —> M and F v : P h —* M by 

B v (u fc ; v h ) := By(u h -v h ) + /3P(u h )P(v h ), (31) 

F v (y h ) := Fy(vh) + /?(1, (1/2 — JC)u 0 )yP(v h ). (32) 

Then the stabilized FVM-BEM coupling reads: find G Ph, such that 

By (u h ; v h ) = F v (y h ) for all w h G 74. (33) 


Remark 11. The discretized version of the stabilized FEM-BEM coupling reads with the 
stabilized weak form (fTTD : find u^fem £ Ph such that 


B (u h ,FEM) v h ) = F(v h ) for all v h G P h . 

See also Remark [?| 

In the spirit of Lemma |3] and |AFF + 13l Theorem 14] we can state the equivalence of the two 
presented FVM-BEM formulations. 


Lemma 12. The FVM-BEM coupling (127|) and its stabilization (j33j) are equivalent. The 
statement is also true if we replace Ay by Ay in the corresponding bilinear forms. 

Proof. In case of fd — 0 the two formulations are obviously the same. Thus we only have to 
consider j3 — 1. If = (v,h, 4>h) is a solution of (127)) . testing with = (0,1) it follows that 

P{ U/») = (1, (1/2 - K)u h + V(f h ) r = (1, (1/2 - JC)u 0 ) r , (34) 

which means that we can add the stabilization term to (1271) to get the stabilized version (135)) . 
Reciprocally, testing (1331) with v h = (0,1), it follows that 

P(u,)(l + (1, Vl) r ) = (1, (1/2 - /C)n 0 )r(l + (1, Vl) r ). 

Since the single layer operator is coercive, (133)) follows and we can eliminate the /5-dependent 
term in (1331) to get (1271) . Note that the proof is independent of the particular choice of the 
finite volume bilinear form, and it therefore holds for Ay as well. □ 

The idea of our analysis is to estimate the difference of the stabilized FEM-BEM coupling 
and the stabilized FVM-BEM coupling. For that we need the following two estimates, which 
are standard in the context of FVM |ELL021. ICha02j with the above constructed dual mesh, 
but here extended to the coupling problem. 


Lemma 13. For the difference of the right-hand side of and (133)) . there holds 

I F(v h ) - Fyfv h )\ < c( Y, ^||/|U 2 (T)||V^|| L 2 (t) 

TeT 

+ X/ hE*\\to ^ io||L2(B)||Vt;/ l ||L2( T£ ,)^ 

EeEr 


(35) 


for all w h = (vhAh) £ P-h with a constant C > 0, which depends only on the shape regularity 
constant. Here, to is the Ey-piecewise integral mean of t 0 and Te is the element associated 
with E. 
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Proof. It is easy to see that from (J9J) and (|26|) we get 

\F{vh) ~ F v (v h ) | = (f,v h -l£v h )n + (t 0 ,v h -T h v h ) r . 

The Cauchy-Schwarz inequality and (I 20 T) (T 22 |) lead to the assertion. □ 

The next lemma gives us an estimate between the weak and the finite volume bilinear form 
for a function Vh G S l (T). 

Lemma 14. Let us assume that b • n is piecewise constant on Y in , i.e. b • n| pin, e V°(£f n ). 
For all Vh,Wh G S 1 (T) there hold 

| A(vJh,Vh) — Av{Whi Vfi)\ < Cl ^ (jlT\\'Wh\\H 1 (T)\\ r Uh\\H 1 (T) S j , (36) 

Ter 

\A(Wh, Vh) — A^(Wh, Vh)\ < C 2 ^ (hT\\u’h\\H 1 (T)\\Vh\\H 1 (T) S j , (37) 

Ter v 

with constants C\ , C' 2 > 0, depending only on the model data A, b, c, and on the shape 
regularity constant. 

Proof. Let us define := Tfvh G V°(T*). Using integration by parts for A(wh,Vh. ) and 

Av(uih,Vh ) the lines in the proof of |Eral21 Lemma 5.2] show with (120]) 

A(w h , v h ) - A v {w h ,v h ) 

= 22 ((-(div A)\7w h + (div h)wh + b • Vw h + cw h ,v h - v* h ) T 
Ter 

+ ((A - A)Vw h • n, v h - v* h ) E (38) 

E£Et 

- 22 (b-n {w h -w h ),v h -vl) E y 

Ee£ T nr in 

Here, div A is the divergence operator applied to the columns of A, A | E G R dxd is the 
average of A over E, and Wh G V°(£r) is the best L 2 (T) approximation of Wh- With a 
standard approximation argument we prove 

\A(w h ,v h ) - A v (w h ,v h )\ 

< CT(\\Wh\\H 1 (T)\\Vh — v V\\l 2 {T) + 22 ^E\\^W h \\ L 2^\\v h — V^\\ L 2^ 

Ter E&£t 

+ "22 W Wh ~ ™h\\L 2 (E)\\Vh ~ V*h\\L 2 { E 2) 

Ee£ T n£p 

with a constant C E > 0 , which depends only on the shape regularity constant, and the model 
data A, b, and c. The standard scaling inequalities ||< Ch~2^ 2 \\ Vw/ 1 ||l 2 (t) and 

||Wh — Wh\\i 2 { E ) < Ch^ 2 ||Vw/j||i 2 ( T ), together with and (]2l]) - (]22|) . prove f]36]) . To prove (1371) 
we write 

I A(w h ,v h ) - Ay p (w h ,v h )\ < | A(w h ,v h ) -A v (w h ,v h ) \ + \ A v (w h ,v h ) - Ay P (w h , v h )\. 

Note that we can directly apply (136|) for the first and |Eral21 Lemma 6.1] for the second 
difference to show (J37j). □ 
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Remark 15. If A is T -piecewise constant, all parts with A vanish in (1331) because of div A = 
0 and (|2Ql) since Vwh is constant. This is well-known and if b = (0, 0) T and c = 0 there 
even holds A(wh,Vh ) = Ay(wh,vD> see e -S ■ [ BR871 lHac89j . Thus the following analysis also 
holds if A is T-piecewise constant. 

Collecting all the results together we prove: 

Lemma 16. Let us assume that b • n is piecewise constant on T in , i.e. b • n| pin e V°(£f n ). 
For all w h = (wh,£h) £ TLh and v h = (vh,ifh) £ Uh there holds 

\B(w h ;v h ) -B v { w h ;v h ) | < (hr\\w h \\ H i(T) IMI/rRT)) (39) 

TeT 

with a constant C > 0, which depends only on the model data A, b, c, and the shape 
regularity constant. The statement is also true if we replace Ay by Aff in the corresponding 
bilinear forms. 

Proof. We estimate 

|B (w h -v h ) - By (w h ; v h ) I = I A(w h , v h ) - Ay(w h , v h ) - (f h , v h - T h Vh ] r | 

<cj:(h 

TeT 

where we used (1351) and (1201) since fh £ V°(£r)- Using (1371) . the proof with Ay' follows from 
this bound. □ 

Theorem 17 (Stability). There exists H > 0 such that the following statement is valid 
provided that T is sufficiently fine, i.e., h := maxyeT hr < H: Let A m i n (A) > Cyc/4 with the 
contraction constant Cjc € [1/2,1) of the double layer potential. Furthermore, let b • n be 
piecewise constant on T m , i.e. b • n| r i™ £ ’P°(£r n ). Then, there exists a constant CVstab > 0 
such that 

By (v h ; v h ) > CvstabII V ?i ||^ for all v h £ H h . (40) 

The constant CV sta b > 0 depends only on the model data A, b, c, the contraction constant 
Cjc, and the shape regularity constant. The statement also holds if we replace Ay by Aff? in 
the corresponding bilinear forms. 

Proof. From (I50j) we see with C' > 0 

By (v h ; v h ) > B (v h ; v h ) - C"/i||t; h ||^i (n) . 

The stability estimate (071) provides B{v h -,v h ) > C' tab ||vfe||^ with C/ tab > 0, which proves 
the coercivity estimate for h small enough. The proof with Ay’ is the same. □ 

Theorem 18 (A priori convergence estimate). There exists H > 0 such that the following 
statement is valid provided that T is sufficiently fine, i.e., h := m‘ax T( zj- hy < H: Let 
A m i n (A) > Cfc /4 with the contraction constant C ^ £ [1/2,1) of the double layer potential 1C. 
Furthermore, let b • n be piecewise constant on T in , i.e. b • n| pin £ V°(£ “)■ For the solution 
u = (u, (jf) £ PL = R 1 (f2) x H-V 2 (y) of our model problem (00)) and the discrete solution 
u h = ( u h, 4>h) £ Bth = 5 1 (T) x V°(£r)) of our FVM-BEM coupling (133)) there holds 

||u“Uft|| w < C est (h\\f\\ L 2 {n ) + h 1 / 2 \\t 0 -t 0 \\ L 2 (r) + (l + h) v inf ||u - v h || w +/i||u|| w ), 
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where t 0 is the £y-piecewise integral mean of t 0 . The constant C est > 0 depends only on 
the model data A, b, c, the contraction constant Cjc, and the shape regularity constant. In 
particular , if u G H 2 (fl), 0 G H l ^ 2 {£ r), and t 0 G H x ^ 2 {£ r), where 

H 1/2 (£ r ) := {u G L 2 (r) | G H 1/2 (E) for all E G £ r } , 

we have first order convergence 

||u - u h \\ n = 0(h). 

The statement is also true if we replace Ay by Aff in the corresponding bilinear forms. 

In the following proof of Theorem [THJ we write the symbol <, if an estimate holds up to 
a multiplicative constant, which depends only on the model data A, b, c, the contraction 
constant Cjc, and the shape regularity constant. 

Proof. For arbitrary w h = ( v h , if h ) G Ti h we define w h = (w h , ip h ) := u h — \ h G TL h - Then we 
get with (SOD 

||u h - v h \\ 2 Hh ~ &v (u ft ; w h ) - B v (v h ; w h ) 

= F v {w h ) - F(w h ) + B (u; w ft ) - By (v h ; w h ) , 

where we used the finite volume discrete system (1331) and the FEM-BEM bilinear from (TT4|) 
with discrete test functions w h G TLh- Since F v (w h ) — F(w h ) = F v (w h ) — F{w h ) we 
apply (155]) and insert w h to estimate 

Il u h ~ v h \\u ~ ^II/IIl 2 (sT)|| Vwh\\L 2 (n) + hl^Wto — to||L 2 (r)||Vw/ l ||L 2 ( f2 ) 

+ B(u-v h ]w h ) +B(v h ;w h ) - By (v h ; w h ), 

where hs r : = max£ £ f r He- For the second term on the right-hand side we apply the bound¬ 
edness of B and we estimate the last two terms with (139]) . Thus we obtain 

||u/i — v h |& < ^||/||l 2 (£T)|| Vu’/i|| L 2( f2 ) + h^ 2 ||t 0 — ^o||L 2 (r)||Vu;/ l ||i2 ( 'Q) 

+ ||u - v h || w ||w/,|| w + h\\v h \\ H ^n)\\w h \\ H i {n) . 

Finally with ||wh||rri(Q) < ||w ft ||x = ||u h - \ h \\u we get 

i /o _ 

||u ft — Vh\\n ^||/||l 2 (q) + h £r ||t 0 — ^o||L 2 (r) + Il u — v h\\n + 

With ||u fc || H i(n) < ||u - v h \\H + ||u|| w and 

||u - u h \\ H < ||u v h \\ n + \\u h - v h \\ H 

we get the assertion with hs v > h. The proof with Afy is the same. □ 

Remark 19. In |Eral21 see Remark 5.1], where we consider a FVM-BEM coupling with a 
three field coupling approach, we have the constraint <f> G L 2 (T) in the case y(x) = 0 from 
assumption (J 2 ]) to get convergence and an error estimate. Note that this regularity is not 
needed therein to prove existence and uniqueness, see ]Eral2l see Remark 5.2], Furthermore, 
there is also an additional assumption necessary in the case ^(x) = 0 , namely divb + 
c = 0 in and b ■ n = 0 on T in . Thus Theorem 0 which essentially shows existence 

and uniqueness of a discrete solution, and Theorem U3 for our non-symmetric FVM-BEM 

coupling are much stronger than what is available for the three field FVM-BEM coupling. 
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FIGURE 2. The error ||V(w — Uh)\\L 2 {o.) i n the if 1 -semi-norm, the error ||w — 
Uh\\L 2 (n) i n the T 2 -norm, and the conormal error || (f> — <fh\\v i n the V-norm in 
the example in subsection 15.11 for uniform mesh-refinement. 

However , the constraint A min (A) > C^/4 on the eigenvalues of A is not needed for the three 
field FVM-BEM coupling. 


5. Numerical results 


In this section we verify our new coupling with three examples. We stress that in all ex¬ 
periments we consider the discrete FVM-BEM system (I25al) - (125bj) and (1271) . respectively, 
where we replace Ay defined in (1241) by the upwind form Afff defined in (130]) if we use an up¬ 
wind scheme for the convection part. We mention once again, that the equivalent stabilized 
FVM-BEM system (155]) is only needed for theoretical reasons. 

All the numerical experiments are done in MATLAB on a standard laptop with a dual core 2.8 
GHz processor and 16 GB memory. Only the implementation of the matrices resulting from 
the V and /C expressions is done in C using the mex -interface of Matlab [Eral2l lEral3aj . 
As introduced earlier, we use the equivalence of norms \\<f — 4>h\\ 2 H -i/2^ ~ \\4> — 4>h\\v 
(V(4> — 4>h),4> ~ 4>h) r, to calculate the conormal error f — <fh- Then \\<f — f>h\\v leads to 
an approximation of a double integral by quadrature. The details can be found in [EralOl 
IEral21 IEral3a] . In all experiments and in each iteration, T consists of triangles, which are 
up to rotation congruent. In this work we only consider uniform mesh refinement, i.e., we 
divide all triangles by four triangles. 
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FIGURE 3. Interior and exterior solution on an uniformly generated mesh 
with 4096 elements in the example in subsection 15.11 


5.1. Mexican hat problem. We consider the square hi = (—1/4,1/4) 2 . We take the 
exact solution to be u(x i,x 2 ) = (1 — 100x 2 — lOOx^e -50 ^!" 1-2 ^) in the interior domain and 
u e (x i,x 2 ) = log \Jx\ + x\ in the exterior. The diffusion matrix is 

\ ( 10 + cosxi 160xix 2 \ j \ 

160xix 2 10 + sinx 2 J ’ ( ' 41j 

and we take b = (0, 0) T and c = 0. Note that in 17 we have A m ; n (A) = 0.342278 and 
Amax(A) = 10.247271. The right-hand side / and the jumps uq and t 0 are calculated ap¬ 
propriately. We stress that u and u e are smooth in 17 and 12 e , respectively. Therefore, we 
expect a convergence order 0 {h l ) for a first order numerical scheme in the TF-norm, where 
h := maxr g 7 -hy denotes the uniform mesh-size. This corresponds to order 0{N~ 1 ^ 2 ) with 
respect to the number of elements N of 77 The initial mesh T ^ consists of 16 triangles. 
Figure [2] shows the curves of the interior error u — Uh in the 77 1 -semi-norm and L 2 -norm, 
respectively, and the conormal error of </>—</>/, in the V-norm. Both axes are scaled logarithmi¬ 
cally; i.e., a straight line g with slope — p corresponds to a dependence g = 0(N~ P ) = 0(h 2p ). 
The interior H 1 -semi-norm error leads to a convergence order 0(N _ 1 7 2 ), whereas the cor¬ 
responding L 2 -norm error decreases with 0{N~ l ). Thus, the error in 77 1 -norm behaves 
like O (A -1 / 2 ). The convergence of the BEM conormal quantity is optimal in the sense of 
0(1V -3 / 4 ) due to the smooth solution. Altogether we see ||u — u h \\ n = 0(N ~ 1 / 2 ) = 0(h ) 
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Figure 4. The error ||V(w — m/i)||l 2 (o) in the H 1 semi-norm, the error 
||m — Uh\\L 2 (n) in the L 2 norm, and the conormal error \\(p — 4>h\\v in the V 
norm in the example in subsection 15.21 for uniform mesh-refinement. 

with u = (u, (p) and = (uh, <f>h) G Uh, which was shown in Theorem [18] for smooth 
solutions. 

Figure [3] shows the solution in and parts of Fl e . We observe the jump on the coupling 
boundary T and remark that the BEM solution is generated pointwise with the aid of the 
exterior representation formula (J5|) on a uniform grid. For points on the boundary T coming 
from the exterior domain, we use the exterior trace of (1511 . Note that instead of (jTJ) this 
approximated trace reads 

u e ,h\r(x) = —(y4>h){x) + ((/C + |^) (u h - m 0 )) (x) (42) 

for a point evaluation x e T, where ip is the interior angle of the intersection of the two 
tangential vectors in x. 

Remark 20. For this example ”f(x) = 0 from assumption (121) . Thus the analysis needs 
the stabilized bilinear form (l3Tj) with /3 = 1 from (ITTh . In particular, we have the condition 
A min (A) > C/c/4, where Cjc G [1/2,1) is the contraction constant of the double layer potential 
1C. Note that our A with A min (A) = 0.342278 fulfills this constraint. If one replace both 
values of 160 by 165 we get A min (A) = 0.003033 which contradicts the bound. However, the 
experiences (not plotted here) show the right convergence behavior. This confirms similar 
observations for FEM-BEM couplings, e.g. [AFF+13], In particular, the bound seems to be 
a theoretical bound also for our FVM-BEM coupling approach. 
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Figure 5. Interior and exterior solution with a weighted upwinding stabi¬ 
lization on an uniformly generated mesh with 4096 elements in the example in 
subsection 15.21 

5.2. Convection-diffusion problem. We consider the model problem on the square do¬ 
main 12 = (0,1/2) x (0,1/2). We choose a fixed diffusion matrix of A = 0.51, a convection 
field b = (1000xi,0) T and a reaction coefficient c = 0. Note that for this problem we do 
not have an inflow boundary T in and thus (fill) is not needed. For all calculations we use the 
upwind discrete coupling with the weighting function $ defined in (1291) . We prescribe an 
analytical solution 



for the interior domain 12 and 

u e (x i, x 2 ) = log ■sj (aq - 0.25) 2 + ( x 2 - 0.25) 2 

for the exterior domain 12 e . We calculate the right-hand side / and the jumps uo and to 
appropriately. Note that A m ; n (A) = 0.5 and that the problem is highly convection dominated. 
The initial mesh T ^ consists of 16 triangles. In Figure [4] we plot the convergence rate for 
uniform mesh-refinement with respect to the number of elements in T. Since the interior 
and exterior solution are smooth as in the previous example in subsection 15.11 we observe a 
similar convergence behavior, in particular, ||u — = 0 {N~ 1 ^ 2 ) = 0 (h ) with u = (u, (j>) 

and u h = [uh, 4>h) £ 'Hh, which also confirms Theorem HU However, due to the strong 
convection, we have a preasymptotic phase. We want to mention, that without any upwind 
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X\ 


(a) Interior numerical solution without stabiliza¬ 
tion. 



Xl 

(b) Contour lines with full upwind. 


Figure 6 . A convection approximation without upwinding or any other sta¬ 
bilization leads to strong oscillations in 
In 


in the example in subsection 15.31 


(b) we see the transmission effects of the interior and exterior problem 


through a contour line plot. 


stabilization, it is not possible to get a stable solution even for more than 4 million elements, 
which is the last mesh in our calculation. In Figure [5] we plot the interior and exterior 
solution. To resolve the shock at aq = 0.25 better and thus to reduce the effects to the 
exterior domain, one can use adaptive mesh refinement as in |Eral3bj . However, this is 
beyond this work. 

5.3. A more practical example. Our last example is a more practical problem. The model 
can describe the stationary concentration of a chemical dissolved and distributed in different 
fluids, where we have a convection dominated problem in 0 and a diffusion distribution in 
Vt e . Note that the interior is a classical model problem and as described above, the coupling 
with the exterior problem can ‘replace’ the boundary condition, which might be difficult to 
find. Our interior domain 0 = (—1/4, l/4) 2 \([0,1/4] x [—1/4,0]) is the classical L-shape. 
The diffusion matrix A = a I in 0 is piecewise constant and reads 

{ 10 -7 for x 2 < 0 , 

10 -6 for x\ > 0 , 

5 • 10 ” 7 else. 

Additionally, we choose b = (15,10) T and c = 10 -2 . The source is in the lower square, i.e. 
/ = 5 for —0.2 < x\ < —0.1, —0.2 < X 2 < —0.05, and / = 0 elsewhere. We prescribe 
the jumps u 0 = 0 and f 0 = 0. Instead of a logarithmic radiation condition, we impose that 
u = doo + C^l/lxl) and |x| — > oo for an unknown G R. An exterior solution of the Laplace 
equation satisfying this type of asymptotic behavior at infinity must have zero average of 
the normal derivative on T, see [ CS85] , We must add a ^ to the representation formulas for 
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Figure 7. Interior and exterior solution with full upwinding stabilization on 
an uniformly generated mesh with 3072 elements in the example in subsec¬ 
tion 15.31 


the exterior solution (|3]) and fl42l) . respectively, and (UJ) becomes 

u e \r (1/2 —|— Af) , u e |r \^(j) T Uoo* 


Thus we have an additional term (^, a oc )r on the left-hand side of (125b|) and an additional 
equation, which ensures (l,(ph)r = 0 as the counterpart. We use the full upwind scheme, 
i.e. (I28]h for the approximation of the convection term and start with a mesh of 12 triangles. 
This example is similar to the one in |Eral21 Subsection 7.2] but with a smaller diffusion. 
Note that the problem is highly convection dominated and the analytical solution is unknown. 
An interior solution without any stabilization is plotted in Figure 1(a) and shows strong 
oscillations. In Figure £(b) we see the contour lines based on a solution generated on a 


mesh T with 49152 elements. The transport is mainly from the source / / 0 in the left 
lower square in the direction of the convection b. We also can see the interaction with the 
exterior domain, hence, the contour lines are circular. In general, the solution of such a 
problem may have local phenomena such as injection wells. As seen in Figure [7] this leads 
to step layers on the boundary (0,0) to (0,—1/4), due to the convection in this direction 
and the different diffusion coefficient of the interior and exterior problem. Since we consider 
here a domain with a reentrant corner and model data with jumps, it is well known that 
uniform mesh refinement can not guarantee optimal convergence rates, i.e. u / H 2 (Q). An 
adaptive mesh refinement steered through a robust a posteriori estimator could lead to a 
more accurate solution as one can find in a similar example for the FVM-BEM three field 
coupling approach in |Eral3bj . 
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6. Conclusions 


We presented a new FVM-BEM coupling method based on the non-symmetric approach 
to solve a transmission problem, i.e., a convection diffusion reaction problem in an interior 
domain coupled with a diffusion process in an unbounded exterior domain. The resulting 
scheme maintains local flux conservation, also in the case when an upwind scheme for con¬ 
vection dominated problems is used. We showed ellipticity of the continuous and discrete 
system or for some model configurations the ellipticity of their equivalent stabilized system. 
Additionally, we could improve the theoretical elliptic constant from previous works. Note 
that the stabilized FVM-BEM system was only used for theoretical purposes. This allowed 
us to show existence and uniqueness, convergence, and an a priori estimate. We stress 
that for some critical model configurations the assumptions on the data and regularity of 
the unknown solution are weaker than for the comparable three held FVM-BEM coupling. 
Moreover, the non-symmetric approach has less discrete unknowns and thus is computational 
cheaper. Our work gives us a recipe for the coupling of BEM with a non-Galerkin method 
like FVM. Our theoretical results were confirmed by three numerical examples, which illus¬ 
trate the strength of the chosen method in terms of local flux conservation and convection 
dominated problems. 
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